Comprehensive analysis of peroxisome proliferator-activated receptors to predict the drug resistance, immune microenvironment, and prognosis in stomach adenocarcinomas

Background Peroxisome proliferator-activated receptors (PPARs) exert multiple functions in the initiation and progression of stomach adenocarcinomas (STAD). This study analyzed the relationship between PPARs and the immune status, molecular mutations, and drug therapy in STAD. Methods The expression profiles of three PPAR genes (PPARA, PPARD and PPARG) were downloaded from The Cancer Genome Atlas (TCGA) dataset to analyze their expression patterns across pan-cancer. The associations between PPARs and clinicopathologic features, prognosis, tumor microenvironment, genome mutation and drug sensitivity were also explored. Co-expression between two PPAR genes was calculated using Pearson analysis. Regulatory pathways of PPARs were scored using gene set variation analysis (GSVA) package. Quantitative real-time polymerase chain reaction (qRT-PCR), Western blot, Cell Counting Kit-8 (CCK-8) assay and transwell assay were conducted to analyze the expression and function of the PPAR genes in STAD cell lines (AGS and SGC7901 cells). Results PPARA, PPARD and PPARG were more abnormally expressed in STAD samples and cell lines when compared to most of 32 type cancers in TCGA. In STAD, the expression of PPARD was higher in Grade 3+4 and male patients, while that of PPARG was higher in patient with Grade 3+4 and age > 60. Patients in high-PPARA expression group tended to have longer survival time. Co-expression analysis revealed 6 genes significantly correlated with the three PPAR genes in STAD. Single-sample GSEA (ssGSEA) showed that the three PPAR genes were enriched in 23 pathways, including MITOTIC_SPINDLE, MYC_TARGETS_V1, E2F_TARGETS and were closely correlated with immune cells, including NK_cells_resting, T_cells_CD4_memory_resting, and macrophages_M0. Immune checkpoint genes (CD274, SIGLEC15) were abnormally expressed between high-PPAR expression and low-PPAR expression groups. TTN, MUC16, FAT2 and ANK3 genes had a high mutation frequency in both high-PPARA/PPARG and low-PPARA/PPARG expression group. Fourteen and two PPARA/PPARD drugs were identified to be able to effectively treat patients in high-PPARA/PPARG and low-PPARA/PPARG expression groups, respectively. We also found that the chemotherapy drug Vinorelbine was positively correlated with the three PPAR genes, showing the potential of Vinorelbine to serve as a treatment drug for STAD. Furthermore, cell experiments demonstrated that PPARG had higher expression in AGS and SGC7901 cells, and that inhibiting PPARG suppressed the viability, migration and invasion of AGS and SGC7901 cells. Conclusions The current results confirmed that the three PPAR genes (PPARA, PPARD and PPARG) affected STAD development through mediating immune microenvironment and genome mutation.


INTRODUCTION
Gastric cancer (GC) is a malignant digestive system tumor that is most common among patients between 40 and 60 years of age.At present, the mortality of GC is increasing annually, accounting for about 1/4 of all tumor-related deaths worldwide (Fukayama et al., 2020).As a malignant pathological phenotype of GC (Bagaev et al., 2021), stomach adenocarcinoma (STAD) is largely resulted from unhealthy lifestyle, genetic predisposition and Helicobacter pylori infection (Bagheri et al., 2018).Currently, STAD is mainly treated by surgery excision, radiotherapy and chemotherapy (Zeng & Jin, 2022), but the prognosis of STAD patients remains dismal (Hoshi, 2020) due to heterogeneity and metastasis.Though immune and targeted therapies can also improve patients' survival outcomes, we currently face a lack of immune checkpoints and treatment targets specific to STAD (Zeng & Jin, 2022).Therefore, discovering potential prognostic markers to facilitate the development of new drugs for STAD is of great significance (Vyve et al., 2023).
PPARs, a subfamily of the nuclear hormone receptor family, function as ligand-activated transcription factors to regulate various biological processes.Binding of PPARs to agonist and then to the 9-cis retinoic acid X receptor generates a heterodimer, which then combines with specific peroxide proliferative response element to alter the transcription of target genes and exerts corresponding biological effect (Chen et al., 2018).PPAR α, PPAR β/ δ and PPAR γ are the three subtypes of PPARs that differ in tissue distribution, selectivity and sensitivity to ligands.Apart from regulating target gene transcription, the three subtypes of PPARs also regulate various physiological processes including sugar and lipid metabolism and inflammatory response (Berger & Moller, 2002;Monsalve et al., 2013).Lu et al., (2005) found that activation of PPAR γ inhibits the occurrence of STAD in mice through apoptosis induced by pioglitazone and tumor suppressor p53.In gastrointestinal tumor cells, PPAR γ agonists induce the apoptosis of various STAD cell lines and G1 phase cell cycle arrest to suppress the invasion of GC cells and metastasis (Chen et al., 2003;Sato et al., 2000;Takahashi et al., 1999).A low level of PPAR γ in STAD cells enhances fatty acid oxidation (FAO) to promote cancer progression (Ezzeddini et al., 2021), and polymorphism of PPAR γ is closely associated with the peptic ulcer disease (PUD) and Helicobacter pylori infection in STAD (Prasad et al., 2008).Another study reported that the CCL20/CCR6 axis mediates PPAR δ dysregulation to promote STAD carcinogenesis through significantly upregulating

Prognosis analysis
Based on the expression of PPAR family genes quantified by a previous study (Huang et al., 2020), the STAD samples were classified into low-PPARA/PPARD/PPARG and high-PPARA/PPARD/PPARG expression groups using surv_cutpoint function in survminer package (Wang et al., 2020).Differences in overall survival (OS) and progression-free interval (PFI) between two groups were compared using Kaplan-Meier survival curve and log-rank test (Tan et al., 2023) in the R package ''survival'' (Shen et al., 2023).

Co-expression analysis of PPARs and functional enrichment analysis
Co-expression between protein-coding genes and PPARs selected under the cut-off value of |cor|>0.2 and false discovery rate (FDR) <0.05 was calculated by Pearson correlation analysis.Functional annotation of genes (FDR+0.005) to Gene Ontology (GO) items and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was performed using clusterProfiler package (Yu et al., 2012) and the results were visualized using the ggplot2 package (Wickham, 2016).

Tumor immune characteristics associated with the PPARs
In TCGA-STAD dataset, relative abundance of 22 immune cells, such as mast_cells_activated, T_cells_CD4_memory_resting, NK_cells_resting, and macrophages_M0, was calculated by CIBERSORT method.Furthermore, a 29-gene signature representing the main functional components of tumors and some other cell populations was obtained from a previous study (Bagaev et al., 2021).Spearman correlation between different immune factors and PPARs were calculated and visualized as heatmaps using the R package ''heatmap''.The samples were classified into two groups based on the median expression value of the PPAR genes, and the differences of the scores of 22 immune cells between two groups were compared using Wilcox.test.

Genome mutation
Somatic nucleotide variation of TCGA-STAD samples from the Genomic Data Commons (GDC) dataset were called by Mutect2.Tumor mutation burden (TMB), immune cell proportion score (IPS) and Tumor Immune Dysfunction and Exclusion (TIDE) score were calculated by maftools package (Mayakonda et al., 2018).

Drug response prediction
A total of 26 stomach cell lines treated by 175 drugs were acquired from Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/)(Yang et al., 2013).Pearson's correlation analysis was performed to calculate drug sensitivity of cancer cells, and statistically significant PPAR genes were selected under |cor|>0.2and FDR <0.05.AUC of an antitumor drug served as an indicator of drug response.Half of the maximum inhibitory concentration (IC50) value was calculated by pRRophetic package (Geeleher, Cox & Huang, 2014) to reflect drug responses to Docetaxel, Vinorelbine, Paclitaxel, and Cisplatin.

QRT-PCR
Total RNAs of GES-1, AGS and SGC7901 were extracted using TRIzol reagent (Sigma-Aldrich, St. Louis, MO, USA).500 ng of the RNA were used to synthesize cDNA applying the HiScript II SuperMix (Vazyme, Nanjing City, China).QRT-PCR was carried out in ABI 7500 System using the SYBR Green Master Mix (Thermo Fisher Scientific, Waltham, MA, USA).GAPDH was an internal reference.The sequence list of primer pairs for the target genes was as follow:

Cell viability
Cell viability of AGS and SGC7901 was measured by performing Cell Counting Kit-8 assay (Beyotime, Jiangsu, China), according to the protocols.The cells treated with or without si PPARG were separately cultured in 96-well plates at a density of 1×10 3 cells/well.After incubation at 37 • C for 2 h, the OD in each well was measured at 450 nm using a microplate reader.

Transwell assay
Cell migration and invasion were detected by conducting transwell assay.Briefly, the cells (5× 10 4 ) were added to the chambers coated (for invasion) or uncoated (for migration) with Matrigel (BD Biosciences, Franklin Lakes, NJ, USA).The upper layer and lower layer were added with serum-free medium and full DMEM medium, respectively.After incubation, migrating or invading cells were fixed with 4% paraformaldehyde and stained with 0.1% crystalline violet for 24 h.Finally, the number of cells was counted under a light microscope.

Abnormal expression of PPARs across pan-cancer
The expression of PPARA, PPARD and PPARG in a total of 32 types of cancers was determined.We observed that PPARA and PPARG genes were upregulated in most para-cancer samples, while PPARD was upregulated in cancer tissues (Figs.1A-1C).The expression of PPARG was upregulated in STAD.The three PPAR genes were differentially expressed in 32 primary tumor samples (Fig. 1D) and were dysregulated to varying degrees in 26 out of 29 cell lines, including in GC (Fig. 1E).These results demonstrated that the PPAR genes played important roles in cancer progression.

The association between the clinicopathologic features and the expression of the PPAR genes
We analyzed the relationship between the expression of PPARs and clinicopathologic staging in GC.There was no statistical difference in the expression of PPARA, PPARD and PPARG genes in Stage I+II and III+IV (Fig. 2A).PPARD and PPARG had higher expression in G3+G4 than that in G1+G2 (Fig. 2B).In samples with an age >60, PPARG was also upregulated (Fig. 2C).The expression of PPARD was higher in males than in females (Fig. 2D).Moreover, STAD samples were divided into high-expression group and low-expression group according to the median expression value of the three PPAR genes.As shown in Table 1, a high PPARD expression was correlated with males.PPARG was correlated with higher grades and older age.Kaplan-Meier curve demonstrated that

Pathway enrichment analysis
GSVA package was employed to calculate the pathway score based on hallmark gene set from GSEA website for each patient in the TCGA dataset.As shown in Fig. 5A, 23 pathways were enriched to proliferation-related pathways such as MITOTIC_SPINDLE, MYC_TARGETS_V1, E2F_TARGETS, G2M_CHECKPOINT, MYC_TARGETS_V2 in tumor samples and metabolism-related pathways such as XENOBIOTIC_METABOLISM, BILE_ACID_METABOLISM, FATTY_ACID_METABOLISM, HEME_METABOLISM in normal samples (Fig. 5A).These results indicated that STAD patients could benefit from taking cell cycle-related immune checkpoint inhibitors.Moreover, the association between pathways and the expression of the PPAR genes in both tumor and normal samples was determined.In tumor samples, PPARD and PPARG were correlated with EMT-related pathways and proliferation-related pathways, respectively (Fig. 5B).In normal samples, PPARA was negatively correlated with proliferation-related pathways, and PPARG was positively correlated with immune-related pathways (Fig. 5C).Based on these results, it can be reasonably concluded that these PPARG genes contributed to STAD progression through activating different pathways, but the specific mechanism remained to be further studied.

Correlation analysis between TME and the PPAR genes
The relative abundance of 22 immune cells was calculated by CIBERSORT analysis, and then the correlation between the PPAR genes and the infiltration of 22 immune NK_cells_resting and macrophages_M0 (p < 0.05, Fig. 6A).Immune score of 29 TMErelated genes was calculated by ssGSEA.Correlation analysis between the PPAR genes and the TME genes showed that PPARA was significantly negatively correlated with antitumor components, matrix remodeling and tumor proliferation rate, while PPARD was positively correlated with matrix composition (p < 0.05, Fig. 6B).Next, we analyzed the relationship between the PPAR genes and seven immune checkpoint genes, and observed that the three PPAR genes were positively related to SLGLEC15 (p < 0.05, Fig. 6C).Specifically, PPARD and PPARG had the highest correlation with SLGLEC15 (R = 0.21, p = 3.8e 05) (R = 0.32, p = 4.2e 10), while PPARA had the highest correlation with LAG3 (R = 0.15, p = 0.0039) (Fig. 6D).Moreover, the expression of the seven immune checkpoint genes in the high-PPAR expression group and low-PPAR expression group was determined.We found that the level of CD274 and SLGLEC15 was upregulated in high-PPARA and high-PPARD expression groups, while CTLA4 and ALGLECL5 were abnormally expressed in high-PPARG and low-PPARG expression groups (Fig. 6E).Finally, the proportion of 22 immune cells in high-PPAR expression group and low-PPAR expression group was calculated, and only few immune cells such as T_cells_CD8, NK_cells_resting showed differences in their proportion between the two groups (Fig. 6F).

Correlation analysis between genome mutation and the PPAR genes
Analysis on the relationship between genome mutation and the PPAR genes showed that the TMB was positively correlated with PPARA and PPARG (Fig. 7A).At the same time, TMB was higher in the high-PPARA and high-PPARG expression groups than the corresponding low expression groups of the two genes (Fig. 7B).In addition, the immune cell proportion score (IPS) and Tumor Immune Dysfunction and Exclusion (TIDE) score revealed that the patients with high expression of PPARG and PPARA were more likely to benefit from taking immune treatment (Figure S1).TTN, MUC16, FAT2 and ANK3 showed higher mutation frequency in high-PPARA-and high-PPARG expression groups, and DNAH11, HERC2, NIPBL had higher mutation frequency in the high-PPARD expression group (Figs.7C-7E).

Drug sensitivity analysis of PPARs
A total of 26 drug-treated STAD cell lines were obtained from GDSC.Pearson analysis was used to calculate the correlation between drug sensitivity and the expression of the three PPAR genes.We identified 14 PPARA-drug sensitivity pairs (Fig. 8A), and there were 8 drug resistance-related pairs and 2 PPARD-drug sensitivity pairs (Fig. 8B).Specifically, OF-1 targeted BRPF2, OSI-027 targeted MTORC2, LGK974 targeted WNT signaling, EPZ004777 targeted chromatin histone methylation, and GI-6780A targeted metabolism (Fig. 8C).Moreover, analysis on the relationship between IC50 of chemotherapy drugs (Docetaxel, Vinorelbine, Paclitaxel and Cisplatin) and the PPAR genes demonstrated that the PPAR genes were positively associated with Vinorelbine (Fig. 8D).

PPARG promoted the invasion and migration ability of GC cell lines
The results of PCR demonstrated that the mRNA level of PPARG was elevated in the two GC cell lines (AGS and SGC7901) compared to normal gastric epithelial cells GES1 (Figs. 9A, 9B).Consistently, the results of Western blot also showed that the protein level of PPARG was significantly higher in AGS and SGC7901 cell lines (Fig. 9C, 9D).Hence, we hypothesized that PPARG could promote the activity of GC cell lines in vitro.The viability of AGS (Fig. 9E) and SGC7901 (Fig. 9F) was significantly lower after suppressing PPARG expression.Transwell assay was performed to detect the migration and invasion ability of AGS and SGC7901 cells after knocking down PPARG.It was observed that the invasion ability of the two GC cell lines was significantly reduced after PPARG inhibition (Figs.9G-9L).These data indicated that PPARG could significantly promote the viability of GC cell lines.

DISCUSSION
PAR α, PPAR β and PPAR γ are three subtypes of PPARs encoded by independent genes.Previous studies reported that these three genes are closely related to the occurrence and progression of GC.The mRNA and protein expression of PPAR γ is significantly higher in GC than that in normal and paracancer gastric mucosa.Previous studies found that the mRNA and protein expression of PPAR γ silenced by mRNAi in human GC MGC803 cells could inhibit cell proliferation (Ma, Yu & Huai, 2009).Lu et al. (2005) confirmed that  PPAR γ ligand troglitazone could reduce the risk of GC induced by chemical carcinogens in rats.In patient-derived GC cells, PPAR δ knockout significantly suppresses the cancer cell invasion and tumor formation (Song et al., 2020).PPAR α/γ agonist TZD18 inhibits the growth of GC cells by inducing apoptosis (Ma et al., 2019).PPAR α shows a high expression in GC and is negatively correlated with prognosis (Chen et al., 2020b).At present, there is no systematic analysis probing into the specific effects of the 3 PPAR genes on GC.To bridge such a gap, this study analyzed the three genes across pan-cancer including GC using bioinformatics analyses, and found that the expression of the three genes was dysregulated to different degrees.We also observed that GC patients with high-expressed PPARA had better prognostic outcomes, while those with high-expressed PPARG exhibited immune regulation potential to control STAD progress.This suggested that the PPAR genes played important roles in the progression of GC.
Co-expression analysis identified six genes that were closely associated with the three PPAR genes.ASAP2, DNM2 and KIF13B encoding signal transduction proteins play crucial roles in tumor proliferation and invasion.Chen et al. (2020a) showed that high-expressed ASAP2 in GC cell promotes cancer progression through regulating miR-770-5p level, while knockdown of ASAP2 can suppress cisplatin (DDP) resistance in GC treatment (Sun et al., 2021).Overexpressed DNM2 often leads to poor prognosis and enhanced invasive phenotype in many cancers (Trochet & Bitoun, 2021).KIF13B interacts with the tubulin VEGFR2, and KIF13B depletion prevents VEGF-mediated neo-vascularization and endothelial migration (Yamada et al., 2014).EAF2 is a potential tumor suppressor that interacts with the eleven-nineteen lysine-rich leukemia (ELL) protein to inhibit cancer progression of B-cell lymphoma, prostatic cancer and hepatocellular carcinoma (Pascal et al., 2019;Xiao et al., 2008).MFSD9 encodes a transport protein, but its function in cancer has not been reported.Autophagy mediated by transmembrane protein TMEM164 degrades ferritin to promote ferroptosis (Liu et al., 2022;Liu et al., 2023a).These genes were closely related to PPARs and they all have critical functions in cancer progression, indicating the complex roles of PPARs in STAD.
The progression and invasion of GC are closely related to the biological characteristics of GC environment (Rojas et al., 2020).Hypoxic microenvironment, acidosis microenvironment and immune inflammatory response are also involved in the development of GC (Peinado, Lavotshkin & Lyden, 2011;Schlappack, Zimmermann & Hill, 1991;Yang, Meng & Wang, 2021).In our research, the PPAR genes were closely associated with T cell CD4 memory resting, NK cell resting, macrophages M0 and immune checkpoints genes such as PD-1 and PDL1.CD4+ T cell fraction is increased in GC (Lu et al., 2017).PPAR-γ inhibition stimulates the differentiation of M2c macrophage-like (CD206(+) CD163(+) CD16(+)) cells (Zizzo & Cohen, 2015).Another study confirmed that PPAR-induced fatty acid oxidation in T cells can improve the efficacy of anti-PD-1 therapy (Chowdhury et al., 2018), indicating that the expression of PPARs is crucial for immune activation and patients' benefit from immune treatment.
TMB is associated with enhanced clinical immunotherapy response in STAD.Samstein et al. (2019) analyzed the correlation between TMB and OS of cancer patients treated with immune checkpoint inhibitor therapy, and they found a positive relationship between higher somatic TMB (highest 20% in each histology) and better OS in a variety of cancer types.Previous studies also reported that MUC16 mutation is predictive of a better prognosis in GC (Jia et al., 2019;Li et al., 2018).Mutation frequencies of MUC16 and TTN are closely correlated with TMB in GC (Yang et al., 2020).In this study, TMB in STAD was also closely related to PPARA and PPARG and the expression of TTN, MUC16, PPARA, and PPARG was positively related.However, despite interesting findings, the current study still had several limitations.Firstly, all the datasets were public databases and the sample size was relatively small.Secondly, the reliability of our findings should be further validated in vivo and in vitro.In addition, several factors including region differences and the presence of other disease should be considered for further analysis.

CONCLUSION
In summary, this study was the first to comprehensively analyze the expression patterns of PPARA, PPARD and PPARG in multiple cancers including in STAD.The three PPAR genes were closely related to the prognosis, immune microenvironment, and genome mutation of STAD.

Figure 2
Figure 2 The association of the clinical features and PPARs genes.(A) The expression differences of PPARs genes in Stage I+II and Stage III+IV.(B) The expression differences of PPARs genes in Grade 1+2 and Grade 3+4.(C) The expression differences of PPARs genes in Age<60 and Age>60.(D) The expression differences of PPARs genes in female and male.Full-size DOI: 10.7717/peerj.17082/fig-2

Figure 3
Figure 3 The survival analysis between high PPARs genes group and low PPARs genes group.(A) Patients in high PPARA group had a longer overall survival (OS) than that in low PPARA group.(B) There is no statistical difference of PFI between high PPARs genes group and low PPARs genes group.Full-size DOI: 10.7717/peerj.17082/fig-3

Figure 6
Figure 6 The association between Tumor microenvironment and PPARs genes.(A) The correlation analysis between 22 immune cells and PPARs.(B) The correlation of 29 TME gene signatures and PPARs genes.(C) The correlation between immune checkpoint genes and PPARs genes.(D) The most relevant immune checkpoint associated to PPARs genes.(E) The expression levels of seven immune checkpoint genes between high PPARs expression group and low PPARs expression group.(F) The expression levels of 22 immune cells between high PPARs expression group and low PPARs expression group.(*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 and ''ns'' is no significant difference.).Full-size DOI: 10.7717/peerj.17082/fig-6

Figure 7 JiaFigure 8
Figure 7 The association between genome mutation and PPARs.(A) The correlation analysis between TMB and PPARs genes.(B) The difference of TMB in the high PPARs expression group and the low PPARs expression group.(C) Gene mutation difference between the high PPARA expression group and the low PPARA expression group.(D) Gene mutation difference between the high PPARD expression group and the low PPARD expression group.(E) Gene mutation difference between the high PPARG expression group and the low PPARG expression group.(*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 and ''ns'' is no significant difference.).Full-size DOI: 10.7717/peerj.17082/fig-7